home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / DELTAT.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  8KB  |  311 lines
  1. /* DeltaT = Ephemeris Time - Universal Time
  2.  *
  3.  * The tabulated values of deltaT, in hundredths of a second,
  4.  * were taken from The Astronomical Almanac, page K8.  The program
  5.  * adjusts for a value of secular tidal acceleration = -23.8946
  6.  * arcsec per century squared.
  7.  *
  8.  * The tabulated range is 1620.0 through 1994.0.  Bessel's interpolation
  9.  * formula is implemented to obtain fourth order interpolated values at
  10.  * intermediate times.  Outside the tabulated range the program
  11.  * calculates approximate formulae of Stephenson and Morrison
  12.  * or K. M. Borkowski.  These approximations have an estimated
  13.  * error of 15 minutes at 1500 B.C.
  14.  *
  15.  * Input Y is the Julian epoch expressed in Julian years.  Y can be
  16.  * found from the Julian date JD by
  17.  *     Y = 2000.0 + (JD - 2451545.0)/365.25.
  18.  * See AA page B4.
  19.  *
  20.  * Output double deltat(Y) is ET-UT in seconds.
  21.  *
  22.  *
  23.  * References:
  24.  *
  25.  * Stephenson, F. R., and L. V. Morrison, "Long-term changes
  26.  * in the rotation of the Earth: 700 B.C. to A.D. 1980,"
  27.  * Philosophical Transactions of the Royal Society of London
  28.  * Series A 313, 47-70 (1984)
  29.  *
  30.  * Borkowski, K. M., "ELP2000-85 and the Dynamical Time
  31.  * - Universal Time relation," Astronomy and Astrophysics
  32.  * 205, L8-L10 (1988)
  33.  * Borkowski's formula is derived from eclipses going back to 2137 BC
  34.  * and uses lunar position based on tidal coefficient of -23.9 arcsec/cy^2.
  35.  *
  36.  * Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
  37.  * and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
  38.  * Their table agrees with the one here, but the entries are
  39.  * rounded to the nearest whole second.
  40.  *
  41.  * Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
  42.  * Eclipse Maps_, Cambridge U. Press (1986)
  43.  */
  44.  
  45. /* If the following number (read from the file aa.ini)
  46.  * is nonzero, then the program will return it
  47.  * and not calculate anything.
  48.  */
  49. double dtgiven = 0.0;
  50. extern double dtgiven;
  51.  
  52.  
  53. #define DEBUG 0
  54. #define TABSTART 1620.0
  55. #define TABEND 1994.0
  56. #define TABSIZ 375
  57.  
  58. /* Note, Stephenson and Morrison's table starts at the year 1630.
  59.  * The Chapronts' table does not agree with the Almanac prior to 1630.
  60.  * The actual accuracy decreases rapidly prior to 1780.
  61.  */
  62. short dt[TABSIZ] = {
  63. /* 1620.0 thru 1659.0 */
  64. 12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
  65. 8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
  66. 6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
  67. 4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
  68. /* 1660.0 thru 1699.0 */
  69. 3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
  70. 2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
  71. 1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
  72. 1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
  73. /* 1700.0 thru 1739.0 */
  74. 900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
  75. 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
  76. 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
  77. 1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
  78. /* 1740.0 thru 1779.0 */
  79. 1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
  80. 1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
  81. 1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
  82. 1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
  83. /* 1780.0 thru 1799.0 */
  84. 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
  85. 1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
  86. /* 1800.0 thru 1819.0 */
  87. 1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
  88. 1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
  89. /* 1820.0 thru 1859.0 */
  90. 1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
  91. 750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
  92. 570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
  93. 710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
  94. /* 1860.0 thru 1899.0 */
  95. 788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
  96. 161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
  97. -540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
  98. -587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
  99. /* 1900.0 thru 1939.0 */
  100. -272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
  101. 1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
  102. 2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
  103. 2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
  104. /* 1940.0 thru 1979.0 */
  105.  2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
  106.  2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
  107.  3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
  108.  4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
  109. /* 1980.0 thru 1990.0 */
  110.  5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
  111.  5686,
  112. /* Extrapolated values 1991.0 thru 1993.0 */
  113.  5740, 5800, 5860,
  114. /* Wild guess for 1994 */
  115.  5930
  116. };
  117.  
  118.  
  119.  
  120. double deltat(Y)
  121. double Y;
  122. {
  123. double ans, p, B;
  124. int d[6];
  125. int i, iy, k;
  126. double floor();
  127.  
  128. if( dtgiven != 0.0 )
  129.     return( dtgiven );
  130.  
  131. if( Y > TABEND )
  132.     {
  133. /* Morrison, L. V. and F. R. Stephenson, "Sun and Planetary System"
  134.  * vol 96,73 eds. W. Fricke, G. Teleki, Reidel, Dordrecht (1982)
  135.  */
  136.     B = 0.01*(Y-1800.0) - 0.1;
  137.     ans = -15.0 + 32.5*B*B;
  138.     return(ans);
  139.     }
  140.  
  141. if( Y < TABSTART )
  142.     {
  143.     if( Y >= 948.0 )
  144.         {
  145. /* Stephenson and Morrison, stated domain is 948 to 1600:
  146.  * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
  147.  */
  148.         B = 0.01*(Y - 2000.0);
  149.         ans = (23.58 * B + 100.3)*B + 101.6;
  150.         }
  151.     else
  152.         {
  153. /* Borkowski */
  154.         B = 0.01*(Y - 2000.0)  +  3.75;
  155.         ans = 35.0 * B * B  +  40.;
  156.         }
  157.     return(ans);
  158.     }
  159.  
  160. /* Besselian interpolation from tabulated values.
  161.  * See AA page K11.
  162.  */
  163.  
  164. /* Index into the table.
  165.  */
  166. p = floor(Y);
  167. iy = p - TABSTART;
  168. /* Zeroth order estimate is value at start of year
  169.  */
  170. ans = dt[iy];
  171. k = iy + 1;
  172. if( k >= TABSIZ )
  173.     goto done; /* No data, can't go on. */
  174.  
  175. /* The fraction of tabulation interval
  176.  */
  177. p = Y - p;
  178.  
  179. /* First order interpolated value
  180.  */
  181. ans += p*(dt[k] - dt[iy]);
  182. if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
  183.     goto done; /* can't do second differences */
  184.  
  185. /* Make table of first differences
  186.  */
  187. k = iy - 2;
  188. for( i=0; i<5; i++ )
  189.     {
  190.     if( (k < 0) || (k+1 >= TABSIZ) )
  191.         {
  192.         d[i] = 0;
  193.         }
  194.     else
  195.         d[i] = dt[k+1] - dt[k];
  196.     k += 1;
  197.     }
  198.  
  199. /* Compute second differences
  200.  */
  201. for( i=0; i<4; i++ )
  202.     d[i] = d[i+1] - d[i];
  203. B = 0.25*p*(p-1.0);
  204. ans += B*(d[1] + d[2]);
  205. #if DEBUG
  206. printf( "B %.4lf, ans %.4lf\n", B, ans );
  207. #endif
  208. if( iy+2 >= TABSIZ )
  209.     goto done;
  210.  
  211. /* Compute third differences
  212.  */
  213. for( i=0; i<3; i++ )
  214.     d[i] = d[i+1] - d[i];
  215. B = 2.0*B/3.0;
  216. ans += (p-0.5)*B*d[1];
  217. #if DEBUG
  218. printf( "B %.4lf, ans %.4lf\n", B*(p-0.5), ans );
  219. #endif
  220. if( (iy-2 < 0) || (iy+3 > TABSIZ) )
  221.     goto done;
  222.  
  223. /* Compute fourth differences
  224.  */
  225. for( i=0; i<2; i++ )
  226.     d[i] = d[i+1] - d[i];
  227. B = 0.125*B*(p+1.0)*(p-2.0);
  228. ans += B*(d[0] + d[1]);
  229. #if DEBUG
  230. printf( "B %.4lf, ans %.4lf\n", B, ans );
  231. #endif
  232.  
  233. done:
  234. /* Astronomical Almanac table is corrected by adding the expression
  235.  *       -1.915914(Julian centuries from 1955)^2 seconds
  236.  * to entries prior to 1955.
  237.  * Entries after 1955 are referred to atomic time standards and
  238.  * are not affected by errors in Lunar or planetary theory.
  239.  */
  240. ans *= 0.01;
  241. if( Y < 1955.0 )
  242.     {
  243.     B = 0.01 * (1955.0 - Y);
  244.     ans += -1.915914 * B * B;
  245.     }
  246. return( ans );
  247. }
  248.  
  249.  
  250.  
  251.  
  252. /* Routine to fill in time values for TDT and UT.
  253.  * These and the input date, JD are global variables.
  254.  * jdflag is set up on reading the initialization file aa.ini.
  255.  */
  256.  
  257. #include "kep.h"
  258.  
  259. int update()
  260. {
  261. double T;
  262. double deltat();
  263.  
  264. /* Convert Julian date to Julian years re J2000.0.
  265.  */
  266. T = 2000.0 + (JD - J2000)/365.25;
  267.  
  268. switch( jdflag )
  269.     {
  270.     case 0:
  271.         TDT = JD;
  272.         UT = JD;
  273.         break;
  274.     case 1:
  275.         TDT = JD;
  276.         UT = TDT - deltat(T)/86400.0;
  277.         jtocal(UT); /* display the other date */
  278.         break;
  279.     case 2:
  280.         UT = JD;
  281.         TDT = UT + deltat(T)/86400.0;
  282.         jtocal(TDT);
  283.         break;
  284.     }
  285. jtocal(JD);
  286. return(0);
  287. }
  288.  
  289.  
  290.  
  291.  
  292.  
  293. /* Exercise program.
  294.  */
  295. #if DEBUG
  296. main()
  297. {
  298. char s[20];
  299. double ans, y;
  300. double deltat();
  301.  
  302. loop:
  303. printf( "year ? " );
  304. gets(s);
  305. sscanf( s, "%lf", &y );
  306. ans = deltat(y);
  307. printf( "%.4lf\n", ans );
  308. goto loop;
  309. }
  310. #endif
  311.